% a=load('/Users/pmanis/Desktop/06414004.atf', 'ascii')

% define parameters of the test pulse
vc = -10;
tstart = 45.6;
tend = 75.6;

% get and scale the time array
x=a(:,1); % time base
x=x*1000; % scale for sample rate

y=mean(a(:,2:end)');  % get mean of the current trace


xst=find(x>=tstart);
xend=find(x>=tend);
xst = xst(1);
xend = xend(1);
clf;
%plot(x(xst:xend)-x(xst), y(xst:xend), 'ro');
xlim([0, 30]);
i0 = mean(y(1:xst-1));
i1 = mean(y(xst+floor((xend-xst)/2):xend));
rin = (i1-i0)/vc;

Q = sum(y(xst:xend));
[Ipeak, jpeak] = min(y(xst:xend));

xfitstart = xst+jpeak;
xfit = x(xfitstart:xend)-tstart;
plot(xfit, y(xfitstart:xend), 'ro');
a0 = i0;
a1 = i0-i1;
t1a = 2.5;
a2 = a1/2;
t2a = 10;
[FP, CHISQ, NITER] = mrqfit('exponential', ...
    [a0, a1*0.7, t1a, a2, t2a ], xfit, y(xfitstart:xend), ...
    [], [], [-5000, -5000, 0.5, -5000, 2], [5000, 5000, 20, 5000, 50], 400, 1e-6 );
tau1 = FP(3);
tau2 = FP(5);
FP
hold on
yfit = FP(1) + FP(2)*exp(-xfit/tau1) + FP(4)*exp(-xfit/tau2);
plot(xfit, yfit, 'g-');
[a0e, a1e, taue] = expfit(xfit, y(xfitstart:xend))


